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Period ratios in multi-planetary systems discovered by Kepler are 
consistent with planet migration 
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ABSTRACT 

The Kepler planet candidates are an interesting testbed for planet formation scenarios. We 
present results from A^-body simulations of multi-planetary systems that resemble those ob- 
served by Kepler. We add both smooth (Type l/ll) and stochastic migration forces. The ob- 
served period ratio distribution is inconsistent with either of those two scenarios on its own. 

However, applying both stochastic and smooth migration forces to the planets simultane- 
ously results in a period ratio distribution that is similar to the observed one. This is a natural 
scenario if planets form in a turbulent proto-planetary disk where these forces are always 
present. We show how the observed period ratio and eccentricity distribution can constrain 
the relative strength of these forces, a parameter which has been notoriously hard to predict 
for decades. 

We make the source code of our simulations and the initial conditions freely available to 
enable the community to expand this study and include effect other than planetary migration. 
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1 INTRODUCTION 

The number of discovered extrasolar planets is increasing rapidly. 
At the time when this letter was submitted, the number of planets 
has reached 78^ Of those, 275 (35%) are in multi-planetary sys- 
tems with two or more planets orbiting the same star. These systems 
are of particular interest to theorists as they can provide valuable 
information about their formation history. 

The existence of mean motion resonances (MMRs) has been 
confirmed in multiple systems. The most studied planetary system 
in a MMR is Gliese 876. It consists of three gas giants which are 
locked in a tight 1 :2:4 Laplace resonance ( [Rivera et al.|2010| Marcy| 
et al.||200l|>. A large number of studies (e.g. |Lee & Peale||2001| 



2l)02||Sneligrove et al.|200T)|Nelson & Papaloizou|2002| |Beaugel 
& Michtchenko 2003 , Veras 2007 1 suggest that migration driven by 
a variety of mechanisms has played an important role in shaping 
this system. 

This is not surprising as planet migration is a natural outcome 
of the interaction of a planet with the proto-planetary disk that it 
forms in (Goldreich & Tremaine 1980). Both ,/V-body and hydro- 
dynamical models can easily reproduce the observed period ratio, 
the eccentricities and the libration pattern even though the precise 
speed at which migration occurs is still up for debate. Furthermore, 
it is possible to place limits on the strength of additional stochastic 
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forces which might be present in a turbulent proto-planetary disk 
(e.g. |Rein & Papaloizou|2009| l. 

Most of these planets have been discovered with the radial ve- 
locity method. This method is biased towards finding heavy planets 
on close-in orbits. Since 2009, the Kepler spacecraft is monitoring 
over one hundred thousand stars. Kepler has discovered thousands 
of planet candidates which now await confirmation (Batal ha"et al.| 
|2012[ l. These planets have on average a much smaller mass than 
those discovered by radial velocity. Kepler therefore opens a new 
window to test planet formation scenarios. 

In this letter, we apply a model of smooth migration (i.e. 
Type I or Type II depending on disk and planet properties) as well 
as stochastic migration forces to each candidate system with multi- 
ple planets. We show that even a small amount of smooth migration 
produces a distribution of period ratio that is inconsistent with ob- 
servations. At each resonance, planets pile up, leading to distinct 
features in the cumulative distribution function of period ratios. 

When stochastic migration forces are added, these features are 
smeared out. By adding just the right amount, one can retain some 
of the features, leading to a period ratio distribution similar to the 
observed one. Our simulations show a pile-up just outside of the 
exact commensurability. This is also seen in the Kepler data. How- 
ever, it is less apparent in our simulations. This could either be re- 
solved by fine-tuning several parameters or by adding additional 
physics such as tidal damping or gap edges. 
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2 SETUP 



We use the REBOUND code (Rein & Liu|2012) to simulate the or- 
bital evolution of planetary systems. We add dissipative forces to 
the equations of motion which mimic the interactions of a planet 
with a proto-planetary disk (Lee & Peale 2001 1. This setup allows 
us to choose two dimensional parameters for each planet, the mi- 
gration timescale r„ and the eccentricity damping time-scale r e . 
Most Kepler planets are not massive enough to open a gap in a 
minimum mass solar nebula (Weidenschilling 1977, Crida et al. 
|2006||Crida|2009| > and migrate in the Type I regime ( |Ward|1986] ^ 
We therefore choose t„ to be a typical value of the Type I migra- 
tion rate for these planets, 10 3 - 10 4 years. Note that many effects 
such as a non-isothermal disk may change the Type I migration 
rate (|Paardekooper et al. 2010). We set the eccentricity damping 
time-scale to be 10 times shorter than the migration timescale in 
all our simulations, r e = t o /10. This choice has been adopted be 
many authors in the past and has been justified by hydrodynami- 
cal simulations (Kley et al. 2004). Tests have shown that the results 
do not strongly depend on this precise value. We use a 15-th order 
RADAU integrator with a time-step set to 10~ 3 times the innermost 
planet's period. Further tests have shown the the results do not de- 
pend on the choice of integrator or time-step. 

To study the effects of migration, we setup systems that 
closely resemble those observed by Kepler. In fact, we initialize 
the entire set of 364 multi-planetary systems and take the stellar 
mass, the planet periods and the planet radii directly from the pub- 
lished KOI tables. From transit observations only the planet radius 
is known. We follow Fabr ycky et al.| ([2012 ) and assume a simple 
mass-radius power-law to get a reasonably estimate of the planet 
mass: m i>/M e = (>Mb) 2 06 . Each planet is initialized on a circular or- 
bit. All systems are coplanar. This seems a reasonable assumptions 
given that most Kepler planets are expected to be highly aligned 
(e.g. Tremaine & Dong 2012, Sanchis-Ojeda et al. 2012 ). Almost 
all systems are stable for at least 10 4 years when initialized this 
way and we ignore those few (< 1%) that are unstable. We also do 
not take into account that several KOI objects have already been 
confirmed and now have improved orbital fits. We are confident 
that this procedure sets up systems that are indeed similar to the 
real planetary systems, at least in an average sense. We have fur- 
thermore tested that a perturbation of the initial orbital parameters 
(such as starting with a more hierarchical period ratio distribution) 
does not change the outcome of our simulations. 



3 SMOOTH PLANETARY MIGRATION 

We add smooth, dissipative migration forces to the equations of 
motions for the outermost planet in each system. All other planets 
only feel the gravitational forces from the star and the other plan- 
ets. This setup naturally leads to convergent migration. As the outer 
planet moves in, it captures the other planets in resonance. Which 
resonance is chosen depends on the initial position, the migration 
speed and the planet masses (Mustill & Wyatt 2011 1. We stop the 
integrations after 10 4 periods of the outer planet at it's initial loca- 
tion. By that time it has moved in significantly. We tested removing 
the migration forces after a different amount of time but did not see 
any qualitative difference. 

In Figure [T] we show the cumulative distribution of period ra- 
tios of neighboring planets in all observed KOI systems as a solid 



Observed KOI planets 
Migration t a =10 4 years 
Migration i a =10 3 years 




Figure 1. Cumulative distribution of the period ratios in KOI systems with 
multiple planets; solid (red) line: observed period ratios, long dashed line: 
period ratios after migration with T a = 10 4 years, short dashed line: period 
ratios after migration with t„ = 10 3 years. 



(red) curve. Note that there are surprisingly few features near inte- 
ger ratio^] However, by closely inspecting the curve near 3:2 and 
2: 1 period ratios, one can see a deficit of planets exactly in the com- 
mensurability and a slight pile-up just outside. 

The dashed curves show the final period ratios in our in- 
tegrations. We plot the results for two different migration rates, 
r„ = 10 4 years and t = 10 3 years. There are now very clear and 
sharp features that can be associated with resonances. Most plan- 
ets end up in the 3:2 resonance followed by the 2:1 and 4:3 reso- 
nance. This is clearly not consistent with the observed distribution. 
By closely inspecting each of these locations, one can see a tiny 
asymmetry favoring larger period ratios over the exact commensu- 
rability. However, this is nowhere near the very apparent deficit of 
planets in the observed period ratio distribution. 



4 STOCHASTIC PLANETARY MIGRATION 

It is very likely that some kind of stochastic force was acting 
on planets at least during some parts of their past. These forces 
could result form MRI turbulence in the proto-planetary disk (e.g. 
|Rein & Papaloizou 2009, Gress el et al.|201l) . Another possibility 
is the gravitational interaction with a remnant planetesimal disk. 
This scenario also leads to migration containing both a smooth and 
stochastic component ( Ormel et al. 2012 1. Even the interaction with 
other planets can be described as a random walk in certain cases 
(see e.g. Zhou et al. 2007 ; Wu & Lithwick 201 1). Each of these pro- 
cesses is not completely explored and it is hard to estimate the pre- 
cise amplitude (or the diffusion coefficient) of the stochastic forces 
in each of those scenarios. We therefore parameterize the forces in 
this study. 

Stochastic forces are modeled following the procedure de- 
scribed in |Rein & Papaloizou] ( |2009| l. The radial and azimuthal 
component of the stochastic force are modeled independently as 
a Markov process. The forces have an exponentially decaying au- 
tocorrelation function with a finite correlation time. This mimics 



Rein et al. 2012 present a similar plot for radial velocity system. 
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Observed KOI planets 
Migration i a =10 4 years, a=10" 6 
Migration T a =10 4 years, a=10" 5 
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Figure 2. Cumulative distribution of the period ratios in KOI systems with 
multiple planets; solid (red) line: observed period ratios, long dashed line: 
period ratios after migration with stochastic forces with a = 10~ 6 , short 
dashed line: period ratios after migration with stochastic forces with a = 

io- 5 . 

the forces of a turbulent accretion disk. We set the correlation time 
to be half of the orbital period of the outer planet. The strength of 
the forces are measured relative to the gravitational force from the 
central star by the dimensionless parameter a. Rein & Papaloizou 
( 2009 1 estimate the value for a to be ~ 5 • 10~ 6 for small mass plan- 
ets that are embedded in a fully MRI turbulent disk (see their Sec- 
tion 3.1). It is important to point out that this value is not well con- 
strained. In our model each planet is forced stochastically, whereas 
only the outer planet feels a net inward (Type I/II) migration force 
as described in the previous section. 

In Figure [2] we show again the cumulative distribution of pe- 
riod ratios. The solid (red) line is the observed distribution. The 
dashed lines are from simulations where both smooth and stochas- 
tic migration forces are turned on. The long and short dashed lines 
correspond to a = 10~ 6 and a = 10~ 5 respectively. 

One can see that when the stochastic forces are sufficiently 
strong, a = 10~ 5 , every resonant feature in the distribution is 
lost. Furthermore, several planetary systems become unstable and 
several planets get ejected (note that the total number of pairs is 
smaller). The number of ejected planets cannot be easily inferred 
from the observed distribution of (not ejected) planets. Additional 
simulation with no net migration (t„ = <x>, not plotted) show the 
same behaviour. 

For a small level of stochastic forcing, a = 10~ 7 , most reso- 
nant feature remain. This is not shown on the plot as it looks almost 
identical to the long dashed line in Figure [T]This is once again not 
consistent with the observed distribution. 

However, for a more moderate level of stochastic forcing, 
a = 10~ 6 , some of the resonant features remain. Especially the 
features near strong resonances such as 4:3 and 2:1 can still be 
seen. This looks surprisingly similar to the observed distribution. 
There is a pile-up of planets just outside the exact commensurabil- 
ity. However, one should note that this model still produces slightly 
too many planets that are too close to the exact commensurability. 
But the tendency of planets to have slightly larger period ratios is 
reproduced correctly. 

Finally, note that the overall net migration of planets moves 
the distribution to the left. This should not be misunderstood as 



Figure 3. Cumulative distribution of the final eccentricities in runs with 
a = 0, 1(T 6 and 1CT 5 , and r„ = 10 3 and 10 4 . 

a prediction but rather illustrates the fact that we start the plane- 
tary systems at their current location and then migrate them closer 
together. Eventually, one should use an initial distribution that is 
motivated by models of planet growth and then integrated forward 
in time until it resembles the observed distribution. But this kind 
of population synthesis introduces many new parameters and we 
therefore decided not to go down this route. 



5 ECCENTRICITY DISTRIBUTION 

When planets migrate and capture into resonances, their eccentric- 
ities rise. Eventually the excitation from the convergent migration 
and the eccentricity damping from the disk reach an equilibrium 
(Papaloizou 20031. Stochastic forces also cause the planets eccen- 
tricities to grow (Rein & Papaloizou 20091. 

If the combination of the smooth and stochastic migration sce- 
narios, as presented in Section]?] is indeed responsible for the ob- 
served period ratio distribution, then we can make a prediction for 
the eccentricity distribution. The mean eccentricity in our simula- 
tion with parameters t„ = 10 4 years and a = 10~ 6 is {e} ~ 0.01. 
In the simulation with parameters t„ = 10 4 years and a = 10~ 5 we 
have (e) ~ 0.05. 

In Figure[3]we plot the final cumulative eccentricity distribu- 
tion in our simulations. The run with a short migration time-scale, 
t„ = 10 3 years, produces high eccentricity planets. Planets are 
captured into resonance earlier and eccentricities grow faster. The 
simulations without stochastic forcing, illustrated by a short and 
medium dashed curve (blue), show signs of a bimodal distribution. 
This is because planets that do (do not) capture in resonance have a 
high (low) eccentricity. Finally, note that the run with a = 10 5 also 
leads to high eccentricities. This is due to the random excitation and 
not due to resonances. 

Unfortunately eccentricities are not known for most of the Ke- 
pler candidates. Only RV follow-up observations or transit-timing 
variations allow a measurement of the eccentricities. However, the 
eccentricity distribution of multi-planetary systems can be esti- 
mated statistically. Tremaine & Dong (20121 and Moorhea d et al.| 
( |2011| report values of e between ~ 0.1 — 0.25. 

It is unfortunate that the simulation which gives the best fit to 
the period ratio distribution is not the best solution in terms of the 
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Shortest period in system longer than 5 days 
Shortest period in system shorter than 5 days 




Figure 4. Cumulative distribution of the period ratios in KOI systems di- 
vided into two bins: planetary systems where the period of the innermost 
planet is shorter/longer than 5 days. 



mean eccentricity. However, it is worth pointing out that observed 
eccentricities have the tendency to be overestimated and the actual 
eccentricities might be lower. This is a well known effect in radial 
velocity observations (Zakamska et al, ||201 1| | and most likely also 
present in the Kepler sample. For example noise in the data or non- 
transiting planets creating additional transit timing variation might 
lead to an over-estimate of the mean eccentricity. The underlying 
reason for this is that the eccentricity is a positive definite quantity. 
In addition to that, there are multiple possible physical solutions to 
this discrepancy. First, one can fine-tune the eccentricity damping 
and migration timescales. The ratio was fixed in our simulations. 
Further experiments have shown that this ratio indeed has an in- 
fluence on the final eccentricity distribution, but it is not possible 
to push the mean eccentricity all the way up to ~ 0.2. Second, the 
length of our integrations was rather short. Longer integration times 
lead to higher eccentricities due to secular interactions and the still 
active stochastic force. Third, massive, non-transiting planets in the 
systems might also alter the eccentricity distribution. All of these 
effects tend to increase the final eccentricity, thus going in the right 
direction. 

A complete model that takes all these effects into account will 
be able to constrain both smooth and stochastic migration forces. 
This is interesting because Type I migration timescales are hard to 
predict as they depend strongly on local disc properties (see e.g. 
|Paarde kooper et al. 2010). A similar problem exists with stochastic 
forces. Dead zones might lower the amplitudes in some regions of 
the disk. However, a precise modeling of these effects goes beyond 
the scope of this letter. 



6 DISCUSSION 

Our results show that neither smooth nor stochastic planetary mi- 
gration alone can reproduce the observed period ratio distribution 
of multi-planetary systems in the Kepler sample. 

However, a combination of those two effects can create a pe- 
riod ratio distribution which is similar to the observed one. If this 
scenario is true, then we can use the eccentricity distribution of 
Kepler planets to constrain the relative strength of stochastic and 
smooth migration forces. 



Other ideas have been brought up to explain the observed pe- 
riod ratios. Terquem & Papaloizou ( 2007 ) show that migration of 
multiple planets in a disk with an inner edge, together with orbital 
circularization causes strict commensurability to be lost. Similar 
scenarios involving tidal interactions have been studied by Delisle 
|eTaiJ ( [20l2l >; |Batygin & Morbidelli| ( |20T2) ; |Lithwick & Wu[ ( |2012f 
There is one important difference to the migration scenario pre- 
sented here. The inclusion of tides leads a strong radial dependence 
of the effect. Such a dependency is currently not observed in the Ke- 
pler data set. This can be verified in Figure [4] where the cumulative 
distribution is divided into two bins with almost equal number of 
planets. The two bins contain systems where the innermost planet's 
period is shorter/longer than 5 days. If tides were an important fac- 
tor in shaping this distribution, one would expect a bigger effect for 
the bin that contains close-in planets. However, both distributions 
are identical. 

Although this letter focused on the migration scenario, we 
are interested in seeing other ideas being tested in a frame- 
work similar to the one presented in this letter. We there- 
fore decided to make the initial conditions, the integrator and 
all plotting routines freely available. They can be downloaded 
as a tar-file from https://github.com/hannorein/rebound/ 
tree/resonancelocation The relevant files for this project 
are located in the directory problems/resonancelocation. We 
hope that this enables the community to further investigate these 
extremely interesting planetary systems and come up with new and 
maybe even better ideas. 
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